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Abstract.  We  consider  a  model  problem  from  collapse  load  analysis  for  plane  strain.  A 
discretized  version  of  the  kinematic  variational  principle  is  solved  by  means  of  a  Newton- 
type  minimization  algorithm.  The  numerical  solution,  which  represents  a  possible  collapse 
mechanism,  is  computed  for  various  shapes  of  the  domain. 


I.  Introduction. 

This  paper  presents  a  new  computational  method  for  plastic  limit  analysis,  based  on  a 
Newton-type  minimization  algorithm  for  the  kinematic  variational  principle,  whose  solu- 
tion represents  a  possible  collapse  mechanism. 

The  model  of  plastic  hmit  analysis  is  by  now  well  known,  see  e.g.  [3,  7].  The  upper 
and  lower  bound  theorem  [3]  identifies  the  Umit  multipUer  of  a  given  load  as:  (i)  the 
largest  multiple  of  it  which  can  be  equihbrated  by  a  staticcJly  admissible  stress  lying 
within  a  specified  yield  domain,  and  simultaneously  as  (ii)  the  minimiun  of  a  certain 
convex  minimization  problem  among  eJI  velocity  fields  which  do  work  against  the  load  at 
a  unit  rate.  The  functional  analytic  setting  of  these  problems  has  been  clarified  by  Strzing, 
Temam,  ajid  others  [S,  9]. 

A  common  approach  to  the  numericsil  solution  of  limit  analysis  problems  is  to  approxi- 
mate the  yield  domain  by  a  polyhedrzJ  domain  in  stress  space.  Then  the  dual  variational 
principles  become  a  dual  pair  of  linear  prograims,  which  may  be  solved,  for  example,  using 
the  simplex  method  [2].  This  approach  is  not  without  disadvantages,  however:  the  linear 
progreonming  problems  one  arrives  at  are  very  large  etnd  highly  degenerate.  As  a  conse- 
quence, the  development  of  improved  algorithms  for  computing  Umit  multipliers  has  been 
an  area  of  continued  attention,  especiailly  for  problems  of  plane  stress,  plate  theory,  and 
three-dimensional  structures,  see  e.g.  [1,4]. 


The  papers  just  cited  all  use  one  or  another  version  of  linear  prograroxning;  in  particular, 
they  £l1I  require  the  replacement  of  a  curvilinear  yield  domain  by  a  polygonal  approxima- 
tion. A  totally  different  approach  is  suggested,  however  by  the  work  of  M.  Overton  [6].  He 
showed  that  the  variational  problem 


(1.1)  inf       /  |Vu| 


L 

can  be  solved  quite  efficiently  by  mejins  of  a  Newton-type  scheme.  The  success  of  his 
method  is  directly  linked  to  the  fact  that  extremals  of  (1.1)  tend  to  have  Vu  =  0  on  large 
portions  of  Q  [8].  Now  (1.1)  is  precisely  the  kinematic  variational  principle  of  hmit  analysis 
for  a  cylinder  with  cross-section  Q.  which  is  loaded  in  antiplane  shear  by  force  per  imit 
cross-sectional  area  /.  Since  Overton's  scheme  works  so  well  in  this  special  case,  it  is 
natural  to  consider  using  a  similar  method  for  a  more  difficult  plasticity  problem  such  as 
plane  strain. 

This  paper  presents  a  method  analogous  to  Overton's  for  a  problem  of  the  form: 

(1.2)  inf  /  |e(u)|,      n  C  R^ 

subject  to  the  constraint: 


(1.3) 


where  Fi  C  dO.,  u  =  (u],U2),  /  now  takes  values  in  R^,  and  e{u)  =  |(Vti  -|-  Vu')  is 
the  lineair  elastic  strain  associated  with  u.  This  is  the  kinematic  variational  principle  of 
plastic  hmit  anailysis  for  a  two-dimensional  body  of  shape  Cl,  loaded  in  traction  eJong  its 
boundary'  Fi  by  force  per  unit  arclength  /,  (see  figure  1).  The  integral  in  (1.2)  is: 

h{u)=    f  \e{u)\ 
Jn 

(1.4)  =/„(e  '.(")')'. 
which  corresponds  to  the  yield  domain 

(1.5)  B=L  =  {a.,)\Y,<'l<-^ 


In  other  words,  the  dual  of  (1.2-1.3),  which  is  the  associated  static  variationaJ  principle  for 
the  plastic  limit  multiplier,  is  given  by: 

(1.6)  max  {A|3a,  ||a||oo  <  1,^  symmetric  ,  div  cr  =  0  in  n,<7  •  n  =  A  •  /  at  Fi}  , 

where  n  is  the  normaJ  vector  to  fi. 

We  should  note  that  it  is  more  common  to  use  the  Von  Mises  yield  domain: 


in  real  applications,  rather  than  our  choice  B.  The  use  of  an  unbounded  yield  domain 
would  comphcate  the  implementation  of  our  method  considerably:  for  example,  replacing 
B  with  B'  would  have  the  effect  of  adding  the  constraint  div  u  =  0  in  (1.2-1.3).  On  the 
other  hand,  it  would  not  be  difficult  to  adapt  our  method  to  handle  an  anisotropic  yield 
condition,  for  example 

with  a, J  >  0  for  ever>'  t,  j.  In  ajiy  event,  our  (1.2-1.3)  should  be  viewed  as  a  model  plasticity 
problem,  effectively  the  first  test  of  Overton's  approach  in  the  context  of  plane  strain. 
We  consider  several  slightly  different  mechanics  applications;  in  all  of  them  we  have  a 


rectamguleLT  bar  with  two  sj-mmetric  cuts  of  depth  {^)  in  the  y-direction.    The  bar  i 


IS 


loaded  in  tension  along  opposite  sides  parallel  to  the  cuts  (see  figure  1).  In  the  different 
applications  the  shape  of  the  cuts  changes. 

The  Newton  method  works  rather  weD.  In  most  cases  the  solution  has  |e(u)|  =  0  on 
most  of  fi,  in  other  words,  the  collapse  mechanism  consists  of  a  dislocation  along  the  curve 
which  joins  the  two  cuts,  with  the  part  of  17  on  either  side  of  the  dislocation  moving  as 
a  rigid  body.  Though  we  know  of  no  theoreticaJ  reason  for  the  collapse  mechjuiism  to 
have  this  form,  it  is  certsunly  plausible.  No  doubt  the  tendency  to  have  such  a  collapse 
mechanism  is  linked  to  the  success  of  the  algorithm. 

The  paper  will  be  organized  as  follows:  in  the  next  part  the  discretization  of  the  problem 
introduced  in  (1]  (with  infinitesimally  thin  cuts)  will  be  presented.   The  method  and  the 


aJgorithm  are  given  in  part  III.  Part  IV  discusses  the  implementation  of  the  aJgorithm  with 
the  construction  of  the  descent  direction  and  the  method  used  to  test  different  shapes  of 
the  domain  with  the  same  prograon.  In  the  last  part  the  numerical  results  are  given  for 
diiFerent  shapes  of  the  cuts. 


■  1/2 
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figure  1:  the  domain  fi  where  the  cuts  are  infinitesimally  thin. 


II  Discretization. 

Using  the  symmetrj*  of  the  domain  and  forces  the  problem  cam  be  solved  on  the  quadratnt 
0<i<|,0<y<^.We  will  still  call  this  paxt  of  the  domain  fi.  The  constraints  airising 
from  the  symmetry  are: 


(2.1) 


I  "2(1, 


y)      =0     0<y<//2, 
0)     =0     0<x<\. 


We  follow  the  notational  conventions  of  Overton  [6].  Let  N  be  the  number  of  mesh 
points  in  each  direction  and  ^  =  l/N  the  mesh  size.  The  mesh  point  corresponding  to  the 
bottom  of  the  cut  (y  =  1/2)  will  be  indexed  by  /  (see  figiire  2).  The  size  of  the  problem 


(1.2-1.3)  is  2N'^.    The  function  u  in  (1.2-1.3)  is  discretized  using  piecewise  lineax  finite 
elements:  its  value  at  the  {i,jy^  mesh  point  is  (ui(t,i),  122(1,;)). 


u. 


figure  2:  discretization  of  the  domain. 


The  discrete  strain  is  noted  e.jy  where  i,j  =  1,...  A^  -  1  and  p  =  1,2  corresponds  to  the 
lower  or  upper  triangle;  thus  on  a  lower  triangle  we  have  (see  figure  3): 

"i(»  +  l,j)-"i(»,;)  A 


(2.2) 
where 


C.jl  =  ^ 


»^2(»,; +l)-ti2(t,j). 


A  =  T^iMiJ  + 1)  -  «i(t\i)  +  Mi  +  1,;)  -  U2(t', ;■)) 


ij  +  1 


figure  3:  indices  on  a  lower  triangle 


A  change  of  indices  gives  the  discrete  strain  on  the  upper  triangle,  e,j2-  During  the 
computation  only  three  terms  of  the  matrix  of  the  strain  need  to  be  stored;  the  fourth  can 
be  deduced  by  sj-mmetry. 

FVom  (1.4)  and  (2.2)  we  deduce  the  discrete  norm  of  the  strain: 


c.;!  I  = 


(tii(t  +  1,;)  -  tii(t,j))^       (U2(t,j-H)-U2(t,i))^ 


(2.3)  +  r 


^2 

1  (ui(i,;  +  1)  -  ui(tj)  +  ti2(t  +  1,;)  -  ti2(«,;)) 


2  /S2 

the  discretized  version  of  (1.2)  is  thus: 

a 
(2.4)  min  V-|e,jpl,      t',j  =  1, . . .  ,;S^  and  p  =  1,2, 

subject  to  the  constraint  coming  from  the  symmetry  (2.1), 

u,(l,i)     =0     1<J<J/, 


(2.5) 

•  I  "2(1,1)      =0     I  <i<N, 

and  the  constraint  (1.3)  : 

(2.6)  ;9  5^/A-,..,(iV,j)  =  l. 

Although  (2.4),  (2.5)  and  (2.6)  is  a  finite  dimensional  convex  optimization  problem,  it  is 
not  easily  solved  because  the  objective  function  to  be  minimized  is  not  different iable  with 
respect  to  u  if  any  of  the  residuals  |e,jp|  is  equaJ  to  zero.  Note  that  the  residual  vanishes  on 
a  triajigle  when  the  displacement  there  corresponds  to  cm  infinitesimal  rigid  body  motion. 
In  the  numerical  discussion  we  will  call  this  situation  a  "zero-residual." 

A  classiczJ  way  to  solve  a  problem  with  a  non-differentiable  function  is  to  consider  the 
continuously  differentiable  approximation  of  /i,  for  example: 

(2-7)  ^(u(t,j))  =  E^(^?;P  +  ^^)^'^^  >  0- 

Notice  that  h,  is  differentiable  everywhere  provided  that  c^  >  0  (we  shall  actually  let  c^ 
depend  on  the  mesh  size  0).    Unfortunately,  for  a  small  value  of  e^  the  problem  is  very 


ill-conditioned  and  for  a  large  value  of  ep  the  solution  we  reach  using  this  approximation 
is  too  far  from  the  solution  of  (1.2-1.3). 

The  algorithm  presented  below  works  directly  with  /i,  not  with  h^.  However,  it  uses  /»« 
to  find  an  initial  displacement  and  to  attempt  to  verify  optimahty  of  the  finjJ  point,  as 
will  be  explained. 

In  the  following  part  we  will  introduce  the  method  we  use  to  solve  this  problem. 

III.  The  method 

Overton  proposed  in  [5]  an  eJgorithm  for  minimizing  functions  of  the  form  :  sum  of 
the  square  roots  of  sums  of  squares.  The  idea  of  the  algorithm  is  to  restrict  the  objective 
function  to  the  linear  subspace  where  the  zero-residuals  remain  unchanged.  In  this  space, 
the  objective  fvmction  is  continuously  differentiable,  emd  so  it  may  be  minimized  using 
Newton's  method.  Then  of  course  one  must  check  whether  the  "solution"  so  obtjiined  is 
also  optima]  in  the  larger  space.  As  presented  in  [5],  the  algorithm  requires  a  "full  rank 
condition"  which  is  not  vslid  for  the  plastic  Umit  analysis  functional.  However,  Overton 
adapted  the  algorithm  in  [6]  to  solve  (1.1).  He  gave  a  simple  way  to  identify  the  subspace 
where  zero-residuals  remain  unchanged  and  to  construct  the  restricted  objective  function. 

The  algorithm  is  organized  as  follows  :  first,  define  the  subspace  where  zero-residuals 
remain  vmchanged  ;  then  update  the  restricted  objective  function,  its  gradient  and  its 
hessiaxi  ;  and  finally,  check  for  optimahty.  Before  going  into  the  details  of  the  algorithm 
and  its  implementation  we  will  give  some  definitions  and  introduce  the  restricted  objective 
function.  Then  a  way  to  check  the  optimality  will  be  presented. 

III.l.  Definition  of  the  subspace  and  the  restricted  objective  function 

We  will  use  the  notations  of  Overton  in  [6].  In  the  following  the  index  j  will  represent 
(*\jiP)i'iJ  =  1, . . .  ,  A'^  -  1  and  p  =  1,2.  The  objective  function  (1.4)  is  written 

M 

(3.1)  h{xi)  =  Y,0\t,{u)l 

where  u  €  R^'^",  A/  =  2  •  (A'  -  1)^  and  e/u)  =  BjU,whereBj  is  a  (4  x  2N^)  matrix. 
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The  active  set  at  the  point  u  is  the  set  of  triangles  with  zero-residuals  : 

(3.2)  J(u)  =  {;||e»|  =  0}, 

zind  the  active  terms  or  residuals  are  ej(u)  such  that  j  G  J{u).  Let  B  be  the  matrix  whose 
columns  are  the  coefficient  matrix  of  the  active  residuals: 

(3.3)  B=  [B,,,Bj„...  ,Bj,],  where  J(u)  =  {;i,j2,...  ,Jp}, 

and  let  Z  =  Z{u),  the  matrix  whose  columns  span  the  null  space  of  B'.    The  matrix  Z 
satisfies  : 

(3.4)  5'Z  =  0, 
and 

(3.5)  e,(u  +  Zp)  =  e»  =  0,     V;  €  7(u),     Vp  €  R^^' 

Now  for  any  it  we  can  define  the  objective  function  h  restricted  to  the  manifold  defined 
by  uQ)R{Z)  as: 

(3.6)  h,{p,)  =  h{u  +  Zp,). 
Consider  the  gradient  and  the  hessian  of  the  inactive  terms: 

(3.7)  9,{u)  =  V|e,(u)|  =  ^ff^yj  ^  ^("), 


.2,    .   „        B,Bl        B,e,(u)e]iu)B] 


(3.8)  G,(u)  =  V^|e»|  =  ^  -  ^^^-L^Jl^yj  ^  7(.) 

then  the  gradient  and  the  hessian  oi  h^  are  given  by: 


(3.9)  Vh,{p,)  =  Z'g{u  +  Zp,l     where  y=    ^    g„ 

(3.10)  V2;i,(p,)  =  Z'G(u  +  Zp,)Z,     where  G=    ^    G,. 
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111. 2.  Testing  the  optimality  of  the  subspace 

Once  a  displacement  is  found  which  is  optimal  in  the  restricted  space,  we  must  check 
whether  it  is  also  optimal  for  the  entire  (unrestricted)  space.  Because  the  full  rsmk  con- 
dition mentioned  eaxher  does  not  hold,  this  is  difficult.  Following  Overton  [6],  we  only 
check  if  we  can  find  a  lower  solution  in  the  unrestricted  space,  using  the  continuously 
differentiate  approximation  of  h  (2.7).  We  perform  few  iterations  of  this  aJgorithm  with 
the  "solution"  in  the  restricted  space  as  initial  point.  If  a  lower  point  is  found,  we  per- 
form another  optimization  using  different  parameters  to  identify  the  active  set  (see  part 
V  for  details).  Although  this  technique  is  not  guaranteed  to  work,  Overton  [6]  gives  some 
arguments  explaining  why  the  algorithm  is  unhkely  to  converge  to  a  nonoptimal  point. 

111.3.  The  algorithm 

The  main  steps  which  must  be  performed  to  obtain  a  new  iterate  u'"*"'  from  u*  are: 

1.  Choose  an  initial  point  u°  for  which  the  hessian  is  non-singular. 

2.  Identify  the  active  set. 

2.1  Merge  into  the  active  set  sJl  the  triangles  which  have  an  exact  zero- residual. 

2.2  Merge  into  the  active  set  all  the  elements  with  a  very  small  residual  provided  that 
in  doing  so  the  objective  fimction  decreases. 

2.3  Only  if  the  objective  function  does  not  increase  significantly,  merge  into  the  active 
set  all  the  elements  with  small  residuals  and  those  surrounded  by  zero  or  almost 
zero-residual  triangle 

2.4  If  the  merging  in  (3.2)  only  is  effective  go  to  (3.1)  to  check  if  there  are  new  zero- 
residuals  to  merge. 

2.5  Update  Z. 

3.  Perform,  an  optimization  step  in  the  space  u  0  R{Z{u)). 

3.1  Compute  a  descent  direction  using  the  restricted  gradient  and  hessian  with  Newton's 
method  if  the  restricted  hessian  is  well  conditioned,  otherwise  with  the  steepest 
descent  method. 

3.2  Perform  a  line  search  in  that  direction. 
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4.    Termination  tests. 

Stop  if  the  number  of  iterations  is  too  large. 

Check  the  optimality  of  the  subspace  using  the  function  h(  if  the  restricted  gradient 

is  sufficiently  reduced. 

IV  Implementation. 

In  this  peirt  we  will  give  some  details  of  the  implementation  of  the  various  steps  of  the 
algorithm. 

IV. 1.  Choosing  the  initial  point. 

The  problem  is  so  badly  conditioned  that  the  initial  point  must  be  chosen  very  carefully. 
Since  we  only  augment  the  active  set  until  we  find  a  minimiim  in  the  subspace,  the  initial 
point  should  not  have  any  of  its  residuals  equal  to  zero.  To  find  an  initial  point  which 
satisfies  this  condition  we  perform  a  few  iterations  of  an  approxirnate  continuously  differ- 
entiable  problem  using  the  function  /i<  defined  in  (2.7)  with  e^  set  equal  to  10~^;9.  This 
method  has  the  advantage  of  giving  a  point  close  to  a  solution  with  all  residuals  strictly 
positive. 

IV. 2.  Identifying  and  updating  the  active  set. 

When  an  element  is  added  to  the  active  set,  i.e.  when  the  value  of  the  corresponding 
residual  is  small,  the  matrix  Z  and  the  \'alue  of  u  must  be  updated.  We  now  consider  the 
relation  between  elements  of  the  active  set,  in  order  to  be  able  to  update  the  matrix  Z 
efficiently. 

A  zero-residual  element  satisifies: 

|e;(u)|  =  0, 

which  is  equi\'alent,  for  a  lower  triangle,  to: 

'  ^/v/2     (u,(z,;-l)-u,(i-l,;-l))  =  0, 
(4.1)        I   0/V2     (ti2(«-l,j)-U2(«-l,;-l))  =  0, 

.   i9/v^      (til(«  -  1.;)  -  U,(l  -  1,;  -  1)  +  U2{l,J  -  1)  -  U2(:  -  1, J  -  1))  =  0, 
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and  for  an  upper  triangle: 

'^/n/5     (ti,(t,j)-ui(.-l,i))  =  0, 

(4.2)  I   ^/v/2     (ti2(t,j)-U2(t,;-l))  =  0, 

^  ^/v/2     (ui(i,j)  -  ui(i,;  -  1)  +  U2(»',»  -  "2(1  -  !,»)  =  0. 

Here  we  see  the  biggest  difference  between  oux  problem  and  Overton's  :  a  zero-residual 
element  is  determined  by  three  parameters  (an  infinitesimal  rigid  motion  in  R'^)  instead 
of  just  one  as  in  [6].  Let  us  write: 

'  *;i      =ui(ij -l)  =  ui{i-l,j -1), 

(4.3)  I   k2     =ti2(«-l,j)  =  "2(t-l,;-l), 

A     =  tzi(t-  -  1,;)  -  ui(t-  -  1,>  -  1)  =  -ti2(»,i  -  1)  +  "2(«  -  1,;  -  1). 

A  problem  2Lrises  which  has  no  analogue  in  Overton's  algorithm:  how  to  define  the 
paxcimeters  (4.3)  for  a  triamgle  newly  merged  with  em  active  set,  if  this  triangle  meets 
another  triangle  of  the  active  set. 

If  the  two  triangles  share  an  edge  then  the  resolution  of  this  problem  is  easy:  the  three 
parameters  axe  uniquely  defined  by  the  relations  (4.3).  We  will  call  a  zero-residuaJ  zone 
such  a  subset  of  J{u).  It  defines  the  motion  of  a  rigid  part  of  the  domain.  Thus  on  a  zero- 
residual  zone,  u  has  a  constant  value  on  eeich  hne  parallel  to  the  x-axis  and  a  constant 
slope  A  in  the  y  direction.  Symmetrically  xx2  has  a  constant  value  in  the  y  direction  and 
a  slope  —A  in  the  i  direction.  The  displacement  on  the  zero  residual  zone  caxx  be  written 

,.  ..  i  ^li^J)     =      ^1  +  A(j -;o), 

(4.4) 


rt^i(',;)    = 


^2  -  A(i-to), 

where  t],  ^2  are  chosen  such  that  ui{ioJo)  =  ^'i  and  U2(io,io)  =  ^2- 

If,  however,  the  two  triangles  share  only  one  mesh  point  then  the  resolution  of  the 
"problem"  is  not  so  clear.  In  principle  the  slope  A  could  be  different  on  each  zone. However, 
we  expect  that  in  most  cases  (near  the  solution)  the  N-alue  of  A  will  be  the  same  on  all 
the  triamgles  hnked  by  a  mesh  point.  A  physiceJ  explanation  of  this  assertion  is  as  follows 
:  Consider  two  rigid  bodies  fij  and  ^2  linked  by  a  point  Pq  (see  figure  5).  The  energy  is 
null  on  fii  U  U2  although  we  know  from  (4.3)  that  A  can  be  different  on  fij  and  ^2-  But 
then  consider  a  convex  domain  Q  such  that  (fij  U  ^2)  C  n  :  the  energ>'  on  the  domain  Q 
depends  strongly  on  the  value  of  the  two  parameters,  i.e.  if  the  displacement  is  constcmt  on 
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fi\ni  Ufij  but  the  value  of  A  is  different  on  f^i  zuid  Q2,  then  the  global  energy  (J"  |e(u)|dn) 
cannot  be  too  small.  We  thus  make  the  convention  that  if  two  zero-residuzJ  zones  meet  as 
in  figure  4,  then  we  shall  combine  them  into  a  single  zera-residucil  zone  (and  hence  require 
that  u  be  the  same  infinitesimaJ  rigid  motion  on  each).  A  numerical  consequence  is  that 
all  mesh  points  of  the  new  zero-residucil  zone  will  satisfy  the  relations  (4.4)  with  the  same 
value  of  the  pairameters. 


figure  4:  value  of  the  displacement  (u  =  ("1,^2))  on  two 
zero- residual  zones  Hnked  by  a  mesh  point. 


Thus  a  triaxigle  such  as  one  maked  *  in  figure  4  should  satisfy  the  relation  (4.3)  and 
belong  to  the  zero-residuaJ  zone.  As  a  result  of  this  convention  the  globzJ  energy  may 
increase  after  the  merging  of  zero- residual  zones.  In  the  following,  by  a  "forced  merging", 
we  mean  a  merging  which  is  performed  although  the  globed  energy  increases. 

To  insure  the  global  convergence  of  the  algorithm  new  triangles  cire  merged  with  the 
active  set  only  if  the  value  of  the  objective  function  does  not  increase  too  much  as  a 
result,  and  if  we  axe  close  to  a  zero-residucJ,  i.e.  if  their  residucds  are  very  small  or  if  the 
condition  number  of  the  hessian  is  very  \aige.  The  value  of  the  displacement  is  stored  if 
it  corresponds  to  the  lowest  energ>'  ever  computed  by  the  algorithm.    If  the  termination 
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test  is  positive,  i.e.  if  the  gradient  is  suflScently  small,  but  the  objective  function  is  not  the 
lowest  ever  computed,  we  must  perform  another  optimization  from  the  stored  value  of  the 
displacement  with  different  parameters  for  the  merging  operation.  (For  the  det^Lils  see  the 
numerical  part). 


figure  5.  position  of  the  two  rigid  bodies  Qj  and  fij 
in  Q  for  different  \'aJues  of  the  parameters  A. 

IV. 3.  Computing  the  restricted  objective  function. 

We  consider  now  the  representation  of  the  displacement  on  Cl.  Oui  goal  is  to  deduce  the 
form  of  the  matrices  B  and  Z  introduced  in  III.l. 

The  set  of  mesh  points  is  divided  into  lists  such  that  every  element  belongs  to  one  and 
only  one  list.  All  mesh  points  in  Ct  which  belong  to  the  same  zero-residual  zone  will  be 
in  the  same  list  and  a  mesh  point  which  does  not  belong  to  any  zero- residual  will  be 
considered  as  a  list  of  one  element. 

The  displacement  on  a  list  is  determined  by  only  two  parameters  (ki,  Jb2)  if  it  is  a  hst  of 
one  element  and  the  three  parameters  Arj,  k2  and  A  defined  in  (4.3)  if  the  hst  has  more  tham 
one  element.  The  first  element  of  the  Hst,  indexed  by  (t'cjo),  is  such  that  ui(to,jo)  =  ^i, 

i'2(»o,;o)  =  *2- 

The  constraint  arising  from  the  symmetry  (2.1)  can  now  be  expressed  in  the  following 
way:   If  the  displacement  in  the  T-direction  or  the  y-direction  of  one  element  of  a  list  is 
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fixed  to  zero  by  the  constraint  (2.1),  then  the  corresponding  parameter  (Jbi  or  ^2)  will  be 
fixed  to  zero  and  the  given  element  will  be  the  first  of  the  list.  (Note  that  if  there  are  two 
or  more  elements  with  a  fixed  displacement,  the  parameter  A  has  to  be  zero  and  it  does 
not  matter  which  one  is  chosen  as  the  first  one).  If  the  displacement  in  the  i-direction 
and/or  the  y-direction  of  two  or  more  elements  axe  fixed  to  zero,  then  the  value  of  the 
corresponding  parameter(s)  ^i,  ^2  aJ^d  ^  3xe  fixed  to  zero.  If  the  value  of  one  or  more 
parameters  of  a  Hst  is  fijced  to  zero  the  Ust  is  called  a  "fixed  list." 

Using  this  structure  and  the  relation  (4.4)  we  can  deduce  a  simple  form  for  the  matrix 
B  defined  in  (3.3).  The  matrix  B  can  be  divided  into  two  column  blocks  corresponding  to 
uj  and  U2,  and  three  row  blocks  corresponding  to  the  three  parameters  tj,  ^2  8Jid  A. 

The  size  of  the  columns  block  is  fixed  but  the  size  of  the  three  row  blocks  depends  on 
the  number  of  parameters  defined  amd  not  fixed  to  zero.  If  we  call  ni, 11.2,113  respectively 
the  size  of  the  three  row  blocks,  the  size  of  the  matrix  is  ((ni  +  n2  +  na)  x  2N^). 

Taking  into  account  that  the  constraint  (2.6)  must  be  satisfied  in  any  direction  of  search 
p,  we  cein  write  B  in  the  following  way  where  the  first  hst  coresponds  to  the  list  represented 
in  figiire  4  with  A  =  A': 

\. 


0000111121221] 


JB-1.23  t-\.ll    0-12  J  ».2 


.y\^ 


uy/ 


The  laist  column  is  defined  by  the  relation  on  the  list  a: 


X^       f{iJ)MiJ)  =       J2       /(».i)("i('o,Jo)  +  U  -  Jo)A,) 


(4.5) 


=  fci  E       /(^■^')     +  ^'  E      0"  -io)/(i,i) 

\(ij)€  list  <  /  \(ij)e  list  s 


Storing  and  factoring  B  would  be  prohibitively  expensive.  Fortunately,  a  full  rank  matrix 
Z  which  spans  N(B^)  can  be  written  in  the  following  form: 

/ 


V 


«,  s\- 


V 


"T 


e.  0,  0,  V 
^  \ 


\ 


\ 


_  u. 


< 


_  u, 


/ 


^ 


/ 


The  matrix  Z  is  divided  in  blocks  conformed  with  those  of  B'.    The  last  row  of  the 
matrix  insures  that  the  constraint  will  be  respected  by  2uiy  search  direction.    To  insure 
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the  numerical  stability,  the  coefficients  of  the  last  row  defined  in  (4.5),  are  scaled  in  the 
following  way: 

d',  =  -djdt     and     D;  =  -DJdt, 

with  t  such  that       \dt\  =  max  {\di\, .  . .  ,  \dg\}  . 

Notice,  there  Eire  no  rows  corresponding  to  uj  for  the  list  "t";  if  this  Ust  has  more  than 
one  element,  then  the  last  row  will  be  repeated  as  many  times  as  necessary. 

Remark  :  This  technique  can  easily  be  extended  to  the  case  where  Ui  and  u^  axe 
constreiined  over  the  whole  domain  by  adding  the  corresponding  terms  in  the  last  row. 
Only  the  hsts  coid  their  structure  are  stored  ;  thus  the  updating  operation  of  B  and  Z  is 
hmited  to  the  updating  of  the  list. 

IV. 4.  The  optimization  in  the  restricted  space. 

As  for  the  problem  of  minimizing  the  energ\'  in  the  space  u  ©  R{z),  we  first  factor  the 
matrix  Z'GZ  and  estimate  its  condition  number  using  the  Gauss  eHmination  method  as 
implemented  in  LINPACK.  If  the  condition  number  is  not  too  large,  we  then  solve  the 
system  of  equations.  If  the  condition  number  is  greater  than  10^°  we  simply  take  the 
negative  gradient  to  be  the  surch  direction.  We  Ccinnot  use  the  conjugate  gradient  method 
to  solve  the  restricted  Newton  system  because  the  numerical  errors  are  too  large;  thus  we 
must  store  the  whole  matrix  Z'GZ  and  increase  the  number  of  computations  per  iteration. 
Note  that  the  restricted  hessian  is  sparse,  therefore  an  algorithm  that  takes  advantage  of 
this  could  be  used  to  speed  up  the  method. 

The  line  search  is  exactly  the  same  as  in  Overton  [5]. 

IV. 5.   Modifying  the  domain. 

The  structure  of  the  algorithm  actually  makes  it  easy  to  handle  changes  in  the  shape  of 
the  domain  with  only  a  minimal  change  in  the  program.  Consider  a  domain  fi'  C  f^  with 
a  regular  boundary.  We  set  the  energy*  outside  of  Q'  to  zero  by  settine: 

M 
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ro   if  j^Q', 

where     o,  =  < 

'        I  1     if     J  €  fi'. 

In  the  next  part,  we  will  present  numerical  results  for  different  shapes  of  the  domain 
obtained  by  using  this  technique. 

V.  Numerical  results. 

The  fortran  program  used  to  obtain  the  results  was  written  by  modifying  and  extending 
the  program  of  Overton  used  to  obtain  the  results  reported  in  [6].  It  was  run  on  a  Vax 
11/780  at  the  Courant  Institute  Mathematics  emd  Computing  Laboratory.  Double  preci- 
sion arithmetic,  i.e.  approximately  16  decimal  digits  of  accuracy,  was  used  throughout. 

We  tried  several  problems  including  the  specific  one  above.  We  considered  different 
shapes  of  the  cuts:  infinitesimal  by  thin  cuts  £is  in  [1]  (figure  6)  ;  large  rectangular  cuts, 
two  mesh  sizes  wide  (figure  7)  ;  and  triangular  cuts  (figure  8).  We  studied  each  of  these 
problems  for  three  different  lengths  of  the  cuts:  /  =  3»  9'  I " 

To  keep  the  symmetry  of  the  continuous  problem  we  have  used  a  symmetric  discretiza- 
tion. Thus  for  the  first  and  second  shape  we  have  solved  the  discrete  problems  for  two 
different  discretizations  (see  figure  9).  In  the  third  problem  only  one  discretization  is 
possible.  In  all  the  tests  presented  here,  we  consider  only  the  constant  loetd  (1.3): 

(6.1)  nx,-^->       '     '  =  5'-3^!'2'- 


^■^'={      0 


0     elsewhere. 

Computational  experience  indicates  that  the  specific  solution  we  reach  and  the  speed  of 
the  method  depends  strongly  on  the  value  of  the  parameters  used  in  step  2  of  the  algorithm: 
f.MCH  which  determines  when  a  residual  is  considered  exactly  active  (i.e.  to  be  a  zero- 
residual).  fACT  which  determines  when  an  attempt  is  mzide  to  meJce  an  almost  active 
residueJ  exactly  active,  ecoND  which  determines  if  the  condition  number  is  sufficiently 
large  to  authorize  a  "forced  merging"  defined  in  rV.2,  and  epoR  which  is  used  when  the 
merging  is  forced  (the  smaller  residual,  while  non-zero,  must  be  smaller  than  epoR  *cact)- 

If  the  termination  test  is  successful  but  the  value  of  the  objective  function  is  not  the 
smaJlest  computed,  the  algorithm  starts  again  with  smaller  \-alues  for  the  psu-ameters  fpoR, 
e.ACT,  CMCH  and  a  larger  \-alue  for  the  parameter  «coND,  with  the  computed  "solution"  as 
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an  initiad  point.  Note  that  a  merging  can  be  "forced",  i.e.  the  objective  function  increased 
after  the  merging,  although  the  residual  was  "exactly  active",  i.e.  |e,l  <  cmch- 
The  paranaeters  are  updated  in  the  following  way: 

«MCH  =      cmch  X  10"  , 

^ACT  =       fACT  ^  10~    , 
fFOR  =        fPOR  X  10"    , 
CCOND  =        fCOND  X  10. 

The  usual  value  of  the  condition  number  is  10®;  a  "large"  value  is  typically  10^^.  If  it  is 
larger  than  10'°  a  steepest  descent  step  is  performed  instead  of  a  Newton  step. 

The  algorithm  is  not  very  sensitive  to  the  value  of  €pg,  which  determines  if  the  gradient 
is  sufficiently  reduced;  it  has  been  fixed  to  10~®.  The  parameters  cline  and  ETA  used  in 
the  line  seztrch  are  the  same  as  in  Overton's  problem  [5,  6]. 

The  tables  I  and  II  summarize  the  results  for  different  problems  and  two  different  initial 
values  of  the  parameters  fACTi^MCH  and  epc-  In  the  last  column  of  the  tables  is  given  an 
estimation  of  the  speed  of  the  algorithm: 


Q 

E 
1=1 


mi 


3 


(5.2)  EST       =      ^    3    , 

Qml 

where  mi  is  the  size  of  the  restricted  space  at  the  tth  iteration,  mo  is  the  size  of  the  initial 
problem  and  Q  the  number  of  iterations.  This  estimator  does  not  take  into  account  the 
number  of  iterations  needed  and  the  iteration  where  the  steepest  descent  method  is  used. 
In  aJl  tests  performed  the  number  of  mesh  points  N  is  12. 

Physical  intuition  suggests  that  the  solution  in  all  cases  should  be: 


(5.3)  ui(3-,y)  =  < 


0  x  =  0,     -|<y<i, 

[  -"1         -^<x<0,      -i<y<i. 


^      .  1111 

u.(x,y)  =  0      -2<-<2'      -2^^^2- 
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The  solutions  we  found  are  very  close  to  this  intuitive  solution.  Note  that  the  form  of 
the  solution  is  part  of  the  success  of  the  eJgorithm  ;  aimong  other  solutions,  this  algorithm 
alw'ays  seems  to  choose  this  type  of  highly  discontinuous  solution. 

We  represent  the  solutions  (figures  10,  11  and  12)  with  the  same  method  as  used  by 
Christiansen  in  [1]:  the  position  {di,d2)  of  a  mesh  point  for  the  optima]  displacement  is 
obtained  by: 

<ili^x,yj)=         Z,  +QUi(l,,t/j) 

<^2(ii,yj)=     Vi +  ou2(zi,yj) 

Our  graphics  axe  obtained  by  using  q  =  5/6.  In  figure  10  the  slight  difference  from  the 
intuitive  solution  is  mainly  due  to  the  discretization.  In  figure  13  we  represent  the  value 
of  the  displacement  after  13  iterations  of  a  Newton  method  applied  to  the  approximated 
problem  (2.7),  using  shape  1,  discretisation  1  with  the  peu-ameter  values  /  =  |iC  =  0.1/9. 
The  vcdue  of  the  projected  gradient  is  3.9  •  10~'^  and  the  vadue  of  the  objective  f.unction 
is  0.944.  Note  that  the  veJuc  of  the  objective  fimction  of  the  same  problem  in  table  I  is 
0.719. 

Conclusion. 

In  this  paper  we  have  described  a  method  to  solve  the  problem  (2.4),  a  discretization 
of  (1.2-1.3).  Although  the  problem  is  more  complex  than  M.  Overton's  problem  [6],  the 
method  hajidles  the  ill-conditioning  of  the  hessiein  very  well,  and  allows  one  rather  quickly 
to  find  a  better  solution  tham  that  obtained  using  a  continuously  different iable  approxima- 
tion of  the  objective  function  (2.7).  The  success  of  the  aJgorithm  is  presumably  linked  to 
the  structure  of  the  solution  (5.3),  which,  though  highly  intuitive,  has  not  yet  been  proved 
(except  numerically)  to  be  optimal. 
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figure  10:   deformation  of  the  domain  (shape  1) 


initial  structure  optimal  deformation  (1=1/2) 

figure  11;   deformation  of  the  domain  (shape  2) 
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initial  structure  optimal  deformation  (1=1/2) 

figure  12:   deformation  of  the  domain  (shaoe  3) 


figure  13  :  value  of  the  displacement  after  13  iterations 

of  the  approximate  problem. 

(shape  1,  1  =  1/2, £=0.1fl) 

value  of  the  objective  function  :  O-^^^^?., 
value  of  the  projected  qradient  i  3.9  10 


Table  II:  Numerical  results  for  the  different  problems. 
iMCH  =  10-^^  ey,cT  =  5  X  10-^Pn,  (PG  =  10-'^  ^  =  12 


Discre- 

Shape 

/ 

Final 

Final  Value 

No.  of 

Iterations 

EST    1 

tization 

o/Q' 

value  of  f 

of  norm  of 
the  gradient 

Newton 

Steepest  Method 

f 

1 

1/3 

.462064 

2.9  X  10-^^ 

26 

26 

17% 

2 

1/3 

.495431 

7.4  X  10-^« 

21 

5 

28% 

1 

1/2 

.7191869 

2.7  X  10-" 

36 

1 

23% 

2 

1/2 

.752561 

0 

26 

1 

35% 

1 

2/3 

.9751832 

6.9  X  10-^* 

28 

9 

55% 

2 

2/3 

1.009578 

1.3  X  10-^ 

32 

5 

35% 

1 

2 

1/3 

.385694 

0 

22 

84 

13% 

2 

2 

1/3 

.385694 

1.3  X  10-^* 

19 

49 

26% 

1 

2 

1/2 

.642824 

0 

26 

3 

31% 

2 

2 

1/2 

.642824 

5.2  X  10-'* 

23 

4 

36% 

1 

2 

2/3 

.899287 

5.7  X  lO-''' 

45 

5 

38% 

2 

2 

2/3 

.899954 

0 

21 

1 

51% 

1 

3 

1/3 

.899954 

0 

28 

0 

27%. 

1 

3 

1/2 

1.028518 

0 

31 

0 

26% 

1 

3 

2/3 

1.157083 

0 

30 

1 

45% 

(*;  the  optimization  was  interrupted  because  the  line  search  could  not  find  a  lower  point.) 
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